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ABSTRACT 

Aperture array interferometers, such as that proposed for the Square Kilometre Array (SKA), 
will see the entire sky, hence the standard approach to simulating visibilities will not be appli- 
cable since it relies on a tangent plane approximation that is valid only for small fields of view. 
We derive interferometric formulations in real, spherical harmonic and wavelet space that in- 
clude contributions over the entire sky and do not rely on any tangent plane approximations. 
A fast wavelet method is developed to simulate the visibilities observed by an interferom- 
eter in the full- sky setting. Computing visibilities using the fast wavelet method adapts to 
the sparse representation of the primary beam and sky intensity in the wavelet basis. Conse- 
quently, the fast wavelet method exhibits superior computational complexity to the real and 
spherical harmonic space methods and may be performed at substantially lower computa- 
tional cost, while introducing only negligible error to simulated visibilities. Low-resolution 
interferometric observations are simulated using all of the methods to compare their perfor- 
mance, demonstrating that the fast wavelet method is approximately three times faster that the 
other methods for these low-resolution simulations. The computational burden of the real and 
spherical harmonic space methods renders these techniques computationally infeasible for 
higher resolution simulations. High-resolution interferometric observations are simulated us- 
ing the fast wavelet method only, demonstrating and validating the application of this method 
to realistic simulations. The fast wavelet method is estimated to provide a greater than ten-fold 
reduction in execution time compared to the other methods for these high-resolution simula- 
tions. 

Key words: techniques: interferometric - methods: numerical - cosmology: observations. 



1 INTRODUCTION 

The next generation of interferometers, such as the Square Kilome- 
tre Array (SKA), will have to overcome a number of challenging 
imaging issues. Key among these is the question of how to deal 
with very large fields of view, both in terms of the forward problem 
of simulating observed visibilities and also in terms of the reverse 
problem of reconstructing wide field images. The reverse wide 
field imaging problem has been tackled traditionally by faceting 
the sky into a number of regions which are sufiiciently small that 
the standard tangent plane approximation to Fourier imaging is 
possible (Cornwell & Perley 1992; Greisen 2002). More recently, 
Corn well et al. (2005) have introduced the w-projection algorithm, 
providing an order of magnitude speed improvement over facet- 
based approaches. The forward wide field imaging problem is an 
issue which arises when observing with interferometric aperture 
arrays, such as those proposed for the low frequency instrument of 
the SKA, and which has yet to be resolved. 

The final configuration and system design of the SKA are still 
under active development. The final design will be dependent on 
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the result of simulations. Such simulations are created specifically 
to assess different array configurations, whilst taking into account 
a range of sky models and possible error contributions. These sim- 
ulations are not only important in determining the response of the 
interferometer to the real sky but also in establishing the dynamic 
range of an observation when faced with bright sources in the side- 
lobes of the primary beam. As a consequence, simulating the re- 
sponse of an aperture array interferometer correctly is as vital to the 
design studies of such instruments as it is important for devising an 
effective method of reducing the data when they finally arrive. 

Such simulations are relatively simple for interferometers 
with small fields of view. Small fields can be well approximated 
as planes tangent to the celestial sphere. This approximation al- 
lows the visibilities observed by an interferometer to be related to 
the tangent plane image through the Fourier transform. It is then 
straightforward to simulate visibilities and reconstruct images. In 
the case of aperture arrays, which may see the entire hemisphere, 
the operation is not so trivial. A tangent plane approximation is 
obviously inappropriate for such geometries and hence the Fourier 
transform can no longer be used to simulate observed visibilities. 

In this article we relax the small field of view assumption and 
tackle interferometry when considering contributions over the en- 
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tire sky. We derive full- sky interferometry formalisms using a num- 
ber of different representations, including representations in real, 
spherical harmonic and wavelet spaces, and discuss the relative 
merits of each approach. Other authors have considered the spher- 
ical harmonic representation of the interferometry integral relating 
the intensity of the sky to the visibilities observed by an interfer- 
ometer, predominantly in the context of observations of the cos- 
mic microwave background (CMB). White et al. (1999) relate the 
visibilities of observations of the CMB to cosmological quantities 
of interest, such as the angular power spectrum, through contact 
with the spherical harmonics. However, the flat tangent plane ap- 
proximation is still made and hence these results remain restricted 
to small fields of view. Bunn & White (2007) extend the work 
of White et al. (1999) to larger fields of view by using mosaic- 
ing, making the flat- sky approximation for each individual point- 
ing used to construct the mosaic. Ng (2001) was the first to present 
the spherical harmonic representation of the interferometry integral 
while including full-sky contributions, also in the context of com- 
puting the angular power spectrum of observations of the CMB. 
Ng (2001) goes on to discuss implications of this result on proper- 
ties of the angular power spectrum recovered from interferometric 
observations. In a separate piece of work, Ng (2005) discusses the 
recovered full- sky power spectra of CMB experiments with asym- 
metric window functions. Although all of these works do address 
spherical harmonic representations of interferometry, and in some 
cases include contributions over the entire sky, the prevailing con- 
text of these works is implications for the angular power spectrum 
recovered from observations of the CMB. Relatively little atten- 
tion has been paid to full- sky interferometry formulations in the 
context of forward or inverse wide field imaging, which pose im- 
portant problems for next generation interferometers. The purpose 
of this article is to address these issues, focusing particularly on the 
forward wide field imaging problem. 

The remainder of this article is organised as follows. In Sec. 2 
we present the mathematical foundations underlying results derived 
later in this paper, including a discussion of harmonic analysis, 
wavelets and rotations on the sphere. In Sec. 3 we derive representa- 
tions of the interferometric visibility integral in real, spherical har- 
monic and wavelet space, while including full- sky contributions. 
The implications of these results for forward and inverse wide field 
imaging are discussed. In Sec. 4 we present resolution simulations 
of full- sky interferometric observations using all of the methods 
discussed in Sec. 3. Concluding remarks are made in Sec. 5. 



2 MATHEMATICAL PRELIMINARIES 

Before presenting the formulation of interferometry on the full sky, 
it is necessary to outline some mathematical preliminaries. We re- 
view harmonic analysis and wavelets on the two-sphere S^, before 
discussing rotations, which are represented by elements of the ro- 
tation group S0(3). By making all assumptions and definitions ex- 
plicit we hope to avoid any confusion over the conventions adopted. 



2.1 Spherical harmonics 

We consider the space of square integrable functions L^(S^, dQ) on 
the unit two-sphere S^, where &Q.{s) - sin 6 &6 d(p is the usual rota- 
tion invariant measure on the sphere and (6, if) denote the spher- 
ical coordinates of s G S^, with colatitude 6 e [O^n] and lon- 
gitude (f G [0, 2;7r). A square integrable function on the sphere 



F e L^(S^, dQ) may be represented by the spherical harmonic ex- 
pansion 

oo i 
i=0 m=-£ 

where the spherical harmonic coefficients are given by the usual 
projection onto the spherical harmonic basis functions through the 
inner product: 

Fe,n= [ F(5)7*,(5)dQ(5). 

The ^ denotes complex conjugation. We adopt the Condon-Shortley 
phase convention where the normalised spherical harmonics are de- 
fined by (Varshalovich et al. 1989) 



where P'^(x) are the associated Legendre functions. Using this nor- 
malisation the orthogonality of the spherical harmonic functions 
reads 



Y^,n{s) F*^,(S)dQ(s) = S^'Smni' , 



where Sij is the Kronecker delta function. 

We complete this section by noting two identities of which we 
will make subsequent use. Namely, we state the addition theorem 
for spherical harmonics 



,2^+1 



(1) 



and the Jacobi- Anger expansion of a plane wave 

gix.j ^ J^^2£ + 1) i' MM \\y\\) Pi(x • y) , (2) 

where x,j e R^, P^(-) is the Legendre function and j>(-) is the 
spherical Bessel function. 

2.2 Wavelets on the sphere 

Wavelets have proved useful in a wide range of applications due to 
their ability to simultaneously resolve signal content in scale and 
position. In order to perform a wavelet analysis on the sphere, it is 
necessary to extend the ordinary Euclidean wavelet framework to 
a spherical manifold. A number of attempts have been made to ex- 
tend wavelets to the sphere. Discrete second generation wavelets 
on the sphere that are based on a multiresolution analysis have 
been developed (Schroder & Sweldens 1995; Sweldens 1996). Fol- 
lowing the generic lifting scheme proposed in these works, Haar 
wavelets on the sphere for particular pixelisation schemes have also 
been developed (Tenorio et al. 1999; Barreiro et al. 2000). These 
discrete constructions allow for the exact reconstruction of a sig- 
nal from its wavelet coefficients but they may not necessarily lead 
to a stable basis (see Sweldens (1997) and references therein). 
Other authors have focused on continuous wavelet methodolo- 
gies on the sphere (Freeden & Windheuser 1997; Freeden et al. 
1997; Holschneider 1996; Torresani 1995; Dahlke & Maass 1996; 
Antoine & Vandergheynst 1998, 1999; Antoine et al. 2002, 2004; 
Demanet & Vandergheynst 2003; Wiauxetal. 2005; Sanz et al. 
2006; McEwen et al. 2006). Although signals can be recon- 
structed exactly from their wavelet coefficients in these contin- 
uous methodologies in theory, the absence of an infinite range 
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of dilations precludes exact reconstruction in practice. Approxi- 
mate reconstruction formula may be developed by building dis- 
crete wavelet frames that are based on the continuous methodol- 
ogy {e.g. Bogdanova et al. 2005). More recently, filter bank wavelet 
methodologies that are essentially based on a continuous wavelet 
framework have been developed for the axisymmetric (Starck et al. 
2006) and directional (Wiaux et al. 2008) cases. These methodolo- 
gies allow the exact reconstruction of a signal from its wavelet co- 
efficients in theory and in practice. In the full- sky interferometry 
fomalism developed herein, we require a wavelet analysis that al- 
lows the perfect reconstruction of a function on the sphere from 
its wavelet coefficients. Furthermore, we also require an orthogo- 
nal wavelet analysis for reasons that will become clear in Sec. 3. 
The only wavelet methodology on the sphere that satisfies these re- 
quirements is the spherical Haar wavelet (SHW) framework, which 
has the additional advantage of simplicity and is also the most com- 
putationally efficient methodology. We adopt the HEALPix^ pixeli- 
sation of the sphere (Gorski et al. 2005) for the implementation of 
the SHW framework due to its hierarchical nature and ubiquitous 
use in the astrophysical community. 

The description of wavelets on the sphere given here 
is based largely on the generic lifting scheme proposed by 
Schroder & Sweldens (1995) and also on the specific definition 
of Haar wavelets on a HEALPix pixelised sphere proposed by 
Barreiro et al. (2000). However, our discussion and definitions con- 
tain a number of notable diff'erences to those given by Barreiro et al. 
(2000) since we construct an orthonormal Haar basis on the sphere 
and describe this in a multiresolution setting. 

We begin by defining a nested hierarchy of spaces as required 
for a multiresolution analysis (see Daubechies (1992) for a more 
detailed discussion of multiresolution analysis). Firstly, consider 
the approximation space Vj on the sphere S^, which is a sub- 
set of the space of square integrable functions on the sphere, i.e. 
Vj c L^(S^, d^^). One may think of Vj as the space of piecewise 
constant functions on the sphere, where the index j corresponds to 
the size of the piecewise constant regions. As the resolution index 
j increases, the size of the piecewise constant regions shrink, until 
in the limit we recover L^(S^, &Q) as j ^ oo.lf the piecewise con- 
stants regions of S^ are arranged hierarchically as j increases, then 
one can construct the nested hierarchy of approximation spaces 



Vi c ^2 c • • • c Vy c L2(S^ dQ) , 



(3) 



where coarser (finer) approximation spaces correspond to a lower 
(higher) resolution level j. For each space Vj we define a basis with 
basis elements given by the scaling functions (/)j^k ^ Vj, where the 
k index corresponds to a translation on the sphere. Now, let us de- 
fine Wj be the orthogonal complement of Vj in Vj+i. Wj essentially 
provides a space for the representation of the components of a func- 
tion in Vj+i that cannot be represented in Vj, i.e. Vj+i = Vj Wj. 
For each space Wj we define a basis with basis elements given by 
the wavelets ij/jj^ G Wj. The wavelet space Wj encodes the diff'er- 
ence (or details) between two successive approximation spaces Vj 
and Vj+i. By expanding the hierarchy of approximation spaces, the 
highest level (finest) space j - /, can then be represented by the 
lowest level (coarsest) space j - \ and the differences between the 
approximation spaces that are encoded by the wavelet spaces: 



Vj-Vi^^Wj. 



(4) 



Let us now relate the generic description of multiresolution 
spaces given above to the HEALPix pixelisation. The HEALPix 
scheme provides a hierarchical pixelisation of the sphere and hence 
may be used to define the nested hierarchy of approximation spaces 
explicitly. The piecewise constant regions of the function spaces Vj 
discussed above now correspond to the pixels of the HEALPix pix- 
elisation at the resolution associated with Vj. To make the associa- 
tion explicit, let Vj correspond to a HEALPix pixelised sphere with 
resolution parameter A^side - (HEALPix spheres are represented 
by the resolution parameter A^side, which is related to the number of 
pixels in the pixelisation hy N - \2N^\aq^). In the HEALPix scheme, 
each pixel at level j is subdivided into four pixels at level j + 1, and 
the nested hierarchy given by (3) is satisfied. The number of pixels 
associated with each space Vj is given hy Nj - 12 x 4^~\ where 
the area of each pixel is given by Aj = AnlNj - 7t/(3 x 4^~^) (note 
that all pixels in a HEALPix sphere at resolution j have equal area). 
It is also useful to note that the number and area of pixels at one 
level relates to adjacent levels through Nj+i - 4Nj and Aj+i = Aj/4 
respectively. 

We are now in a position to define the scaling functions and 
wavelets explicitly for the Haar basis on the nested hierarchy of 
HEALPix spheres. In this setting the index k corresponds to the po- 
sition of pixels on the sphere, i.e. for Vj we get the range of values 
k = 0, - • • ,Nj - I, and we let Pj^k represent the region of the ^th 
pixel of a HEALPix sphere at resolution j. For the Haar basis, we 
define the scaling function 0^^^ at level j to be constant for pixel k 
and zero elsewhere: 







elsewhere . 



The non-zero value of the scaling function 1 / yA^ is chosen to en- 
sure that the scaling functions (f)j^k for k = 0, - ,Nj - \ do indeed 
define an orthonormal basis for Vj. Before defining the wavelets 
explicitly, we fix some additional notation. Pixel Pj^k at level j is 
subdivided into four pixels at level 7 + 1, which we label Pj+i,kQ, 
Pj+i,ki, Pj+i,k2 and Pj+i,k2^ as illustrated in Fig. 1. An orthonormal 
basis for the wavelet space Wj, the orthogonal complement of Vj, 
is then given by the following wavelets of type m = {0, 1, 2): 



il/\^{s) = [(pj+i,ko(s) + (pj+i,k,(s) - (pj+i,k2(s) - (pj+i,k,(s)]/2 ; 

We require three independent wavelet types to construct a complete 
basis for Wj since the dimension of Vj+i (given by Nj+i) is four 
times larger than the dimension of Vj (the approximation function 
provides the fourth component). The Haar scaling functions and 
wavelets defined here on the sphere are illustrated in Fig. 1. 

Let us check that the scaling functions and wavelets satisfy the 
requirements for an orthonormal multiresolution analysis as out- 
lined previously. We require Wj to be orthogonal to Vj, i.e. we re- 
quire 
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This is always satisfied since for k' k the scaling function and 
wavelet do not overlap and so the integrand is zero always, and for 
k' - kwQ find 

0,,,(s) ilfl,{s) dQ(s) oc J^^ 0- (s) &a{s) = . 
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We also require Wj to be orthogonal to Wy for all j and / . Again, 
if the basis functions do not overlap {i.e. k ^ k') then this require- 
ment is satisfied automatically, and if they do (i.e. k - k') then the 
wavelet at the finer level / > j will always lie within a region of 
the wavelet at level j with constant value, and consequently 

Finally, to ensure that we have constructed an orthonormal wavelet 
basis for Wj, we check the orthogonality of all wavelets at level j: 

where for mi^ m' the positive and negative regions of the integrand 
cancel exactly and for k k' the wavelets do not overlap and so the 
integrand is zero always. Note that in the previous expression the 
final Aj term arises from the area element. The Haar approxima- 
tion and wavelet spaces that we have constructed therefore satisfy 
the requirements of an orthonormal multiresolution analysis on the 
sphere. Although the orthogonal nature of these spaces is impor- 
tant, a diff'erent normalisation could be chosen. It is now possible 
to define the analysis and synthesis of a function on the sphere in 
this Haar wavelet multiresolution framework. 

The decomposition of a function defined on a HEALPix sphere 
at resolution 7, i.e. Fj G Vj, into its wavelet and scaling coefficients 
proceeds as follows. Consider an intermediate level j + \ < J and 
let Fj+i be the approximation of Fy in Vj+i . The scaling coefficients 
at the coarser level j are given by the projection of Fj+i onto the 
scaling functions (pjy. 

Aj,k= [ Fj^,{s)(l>j,k{s)da{s) 

where we call X^k the approximation coefficients since they define 
the approximation function Fj e Vj. At the finest level /, we natu- 
rally associate the function values of Fy with the approximation co- 
efficients of this level. The wavelet coefficients at level j are given 
by the projection of Fj+i onto the wavelets if/Jj^: 

jl,^ Mis) ^l,{s) dais) ^ 

giving 

y\k = {^i+hh + ^i+Ui - ^7+1,^3) 

and 

where we call y^^ the detail (or wavelet) coefficients of type m. 
Starting from the finest level /, we compute the approximation and 
detail coefficients at level 7 - 1 as outlined above. We then repeat 
this procedure to decompose the approximation coefficients at level 
/-I {i.e. the approximation function Fj_i), into approximation and 
detail coefficients at the coarser level J-2. Repeating this procedure 
continually, we recover the multiresolution representation of Fy in 
terms of the coarsest level approximation Fi and all of the detail 
coefficients, as specified by (4) and illustrated in Fig. 2. In general 
it is not necessary to continue the multiresolution decomposition 
down to the coarsest level j = 1; one may choose to stop at the 
intermediate level Jq, where \ Jq < J. 



The function Fj ^ Vj may then be synthesised from its ap- 
proximation and detail coefficients. Due to the orthogonal nature of 
the Haar basis, the approximation coefficients at level 7 + 1 may be 
reconstructed from the weighted expansion of the scaling function 
and wavelets at the coarser level 7, where the weights are given by 
the approximation and detail coefficients respectively. Writing this 
expansion explicitly, the approximation coefficients at level 7 + 1 
are given in terms of the approximation and detail coefficients of 
the coarser level 7: 

^j^iM = ihk + 7% + y],k + ylk)l ; 

= ihk - y\k + y\k - y],k)l ^l^j ; 
^M,k2 = {hk + ylk - y],k - y],k)l ^l^j ; 
^M,k, = ihk - ylk - y\k + y\k)l • 

Repeating this procedure from level 7 = 7o up to 7 = /, one finds 
that the signal Fj gVj may be written as the scaling function and 
wavelet expansion 

Fjis) = £ Aj,,, cPj,,,{s) + X Z Z • (5) 

k=0 j=JQ k=0 m=0 

2.3 Rotations 

Rotations 9^ on the sphere are characterised by the elements of the 
rotation group S0(3), which we parameterise in terms of the three 
Euler angles p = (a,j3,y) e S0(3), where a e [0,2;r), /3 e [0,;r] 
and 7 e [0, In). The rotation of a coordinate vector s by 9^(p) may 
be represented by multiplication of the Cartesian coordinate with 
the 3x3 rotation matrix R(p). We adopt the zyz Euler convention 
corresponding to the rotation of a physical body in a fixed coordi- 
nate system about the z, y and z axes by 7,/? and a respectively, i.e. 
R(p) = R,{a)Ry(fi)R,{y), where R,{§) and Ry{§) are rotation ma- 
trices representing rotations by § about the z and y axis respectively. 
The inverse rotation is given by R~^(p) = Rz(-y)Ry(-/3)Rz(-a). 
We define the rotation of a function F e L^(S^, d^^) on the sphere 
by 

(1^(p)F)Cs) = F(R-'(p)s). (6) 

It is also useful to characterise the rotation of a function on the 
sphere in harmonic space. The Wigner Z)-functions D^^^ip) provide 
the irreducible unitary representation of the rotation group S0(3). 
For our purpose, we merely consider the Wigner Z)-functions to 
represent the rotation of a function on the sphere in harmonic 
space. The rotation of a spherical harmonic basis function may 
be represented by a sum of weighted harmonics of the same £ 
(Brink & Satchler 1999): (n(p)Y„n){s) = ZL-f ^L(P) Y^nis). It is 
then trivial to show that the harmonic coefficients of a rotated func- 
tion are related to the coefficients of the original function by 

inp)F),^ = DL(p) Fen ■ (7) 

n=-£ 

For computational purposes, the Wigner functions may be decom- 
posed as D^^^(a,j3,y) = e~^'"^ ^LO^) where the real polar 
J-functions are defined by Varshalovich et al. (1989). Recursion 
formulae exist to compute the Wigner i/-functions rapidly (Risbo 
1996). 
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Level j + 1 




4>3,k{^) ^j,ki^) ^Iki^) ^Iki^) 



Figure 1. Haar scaling function (pjj^is) and wavelets ^^^(s). Dark shaded regions correspond to negative constant values, light shaded regions correspond to 
positive constant values and unshaded regions correspond to zero. The scaling function and wavelets at level j and position k are non-zero on pixel Pjj^ only. 
Pixel Pj^k at level j is subdivided into four pixels at level j + 1, which we label Pj+i^kQ, Pj+iM , Pj+iM and Pj+i,k3 as defined above. 



Level J - 2 




Figure 2. Haar multiresolution decomposition. Starting at the finest level / (the original sphere), the approximation and detail coefficients at level / - 1 are 
computed. This procedure is repeated to decompose the approximation coefficients at level / - 1 (i.e. the approximation function Fj-\), into approximation 
and detail coefficients at the coarser level J -2. Repeating this procedure continually, one recovers the multiresolution representation of in terms of the 
coarsest level approximation Fj^ and all of the detail coefficients. 
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3 FULL-SKY INTERFEROMETRY 

We formulate full- sky interferometry in this section in real, spher- 
ical harmonic and SHW spaces. Since we treat interferometry in 
the full-sky setting it is necessary to be explicit about the coordi- 
nate systems used and the relations between them. After discussing 
the various coordinate systems that we use, we describe visibility 
computation and image reconstruction for the representation of in- 
terferometry in each space, highlighting the relative merits of each 
representation. These ideas are then extended to incorporate hori- 
zon occlusion and primary beam functions that depend on the in- 
terferometer pointing direction. 

3.1 Coordinate systems 

The complex visibility measured by an interferometer is given by 
the coordinate free definition (Thompson et al. 2001) 

^(m) = ^ A(cr)/(o-)e-^2^"-^ dQ , (8) 

where a - s - Sq\ area element dQ and the vectors cr, s and 
Sq are defined in Fig. 3. We assume that the sky is mapped onto 
the unit celestial sphere, hence all vectors on the sphere are unit 
vectors. The vector u is related to the baseline of the interferometer 
(and is defined explicitly later), the function A defines the primary 
beam of the interferometer and the function / defines the intensity 
of the sky at position cr. In this coordinate free definition of visibil- 
ity, a is essentially the representation of 5 in a coordinate system 
centred on Sq. The translation a - s - sq represents the transfor- 
mation between the global coordinate frame of s and the local co- 
ordinate frame of a. In general, one can transform vectors between 
global coordinates and local coordinates relative to 5o, through a 
rotation by Sq. We now formalise this notion. 

Let us define two coordinate frames. The global coordinate 
frame of the celestial sky is defined by the right handed set of unit 
vectors {wi, W2, n^]. The local coordinate frame relative to Sq is de- 
fined by the right handed set of unit vectors {mi,M2,«3}, where 
aligns with Sq. These two coordinate frames are related by the rota- 
tion = 1^((po, 6o, 0), where (^o, ^o) are the spherical coordinates 
of So in the global coordinate frame. The corresponding 3 X 3 ro- 
tation matrix Rq = R(v^o, 0) transforms a vector between local 
and global coordinates. For example, Sq in global coordinates corre- 
sponds to U3 in local coordinates (by definition) and the two vectors 
are related by Sq = Ro«3- In general, local coordinates are related 
to global coordinates by - Rq^^", where the superscripts 1 and 
n denote local and global coordinates respectively; henceforth all 
vectors and functions contain an 1 or n superscript to denote their 
coordinate frame. The third Euler angle of the rotation mapping lo- 
cal to global coordinates specifies an azimuthal rotation around Sq, 
which is in general arbitrary but fixed. Without loss of generality 
we set this third Euler angle to zero. Fig. 4 illustrates the rotation 
that transforms between local and global coordinates. 

Returning to the visibility function, we may now represent this 
in the coordinate frames defined above. The beam function is most 
naturally represented in local coordinates relative to the pointing 
direction Sq. We denote this function by A\t) e V(S^, dQ). The 
source intensity function is most naturally represented in global co- 
ordinates and is denoted by F(s^) e L^(S^, dQ). We may convert a 
function in global coordinates to a corresponding function in 
local coordinates through the rotation 7^o- 

FXs") = F\Rof) = (^o'^")(^') = F\s') , 
i.e. F^ = In practice, sampled functions on the sphere may 




fl2 



Figure 3. Geometry of an observation of an extended source centred at 
the interferometer pointing direction sq. The area element dQ represents 
the contribution to the visibility integral from point s. The full visibility 
is obtained by summing all such contributions over the sky. Note that the 
sky has been mapped onto the unit celestial sphere, hence s and sq are unit 
vectors. 

be rotated in real space through (6) or alternatively, and more ac- 
curately (since pixelisation artifacts are eliminated), in harmonic 
space through (7). In this coordinate setting the visibility integral 
may be written 

When dealing with the visibility integral it is more convenient to 
represent all functions in the same coordinate frame. In local coor- 
dinates the visibility integral becomes 

^(m)= f A^(5^) c;^o^r)(5V"'^''"'^' di^(^^) 

= f A\t) e-i2^"-^' . (9) 

Here u is the interferometer baseline in local coordinates (an ex- 
plicit expression for u is deferred until the more general formu- 
lation presented in Sec. 3.4). The expression given by (9) is the 
familiar interferometric visibility integral, however in discussing 
this integral in the full-sky setting it has been necessary to make 
the use of diff'erent coordinate systems explicit. To compute visi- 
bilities when including contributions over the entire sky, one could 
simply evaluate (9) by using an appropriate quadrature rule on the 
sphere. The complexity of evaluating this integral directly for a sin- 
gle baseline u scales as 0{N) (recall that N is the number of pixels 
contained in the pixelisation of the sphere). Low-resolution simula- 
tions of full- sky interferometric observations are computed in this 
manner in Sec. 4. 

We also define the coordinate system p - (p, q) on the tangent 
plane of the celestial sphere at the pointing direction. This coordi- 
nate system is used to define the image plane when reconstructing 
images on small patches using the standard Fourier transform ap- 
proach and is obviously a local coordinate frame. It is related to 
the spherical coordinates (^\v^^) of the local [ui, 112,113] frame by 
p - sin cos Lp^ and q = sin 6^ sin ip^. 

3.2 Spherical harmonic space representation 

It is more general and accurate to compute the visibility integral 
in harmonic space since this avoids the need for any quadrature 
rule on the sphere, which would necessarily be pixelisation depen- 
dent and not always exact. Moreover, rotations can be performed 
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Figure 4. Rotation mapping global coordinates of the celestial sky de- 
fined in the frame {h\,h2,h^] to local coordinates defined in the frame 
{mi,M2,M3}, where M3 is aligned with sq. Local coordinates are related to 
global coordinates by = R^^s". 

through (7) in harmonic space more accurately than they can be 
performed in real space. We derive here the spherical harmonic rep- 
resentation of the visibility integral and discuss the implications of 
this representation for computing visibilities and for image recon- 
struction. 

3. 2. 1 Computing visibilities 

Consider the beam-modulated source intensity function 
(A^ • = A\t)l\t). Substituting the spherical harmonic 

expansion of this function into the local coordinate visibility 
integral given by (9), we obtain 

^^{u) = ^ (A^ • 7%^ r e-^2."-^' Yemif) , (10) 

where (A^ • P)^^ are the spherical harmonic coefficients of the 
beam-modulated intensity function. We assume that the beam- 
modulated intensity function is band-limited at ^max so that all 
higher frequency harmonic coefficients are zero, i.e. (A^ • /)^,^ = 0, 
> ^max- In this case sums over the harmonic index i may be 
truncated at ^max- Here, and subsequently, we use the shorthand no- 
tation = 2^^Q Yfm=-£- The integral contained in (10) may be 
evaluated analytically. Noting the addition theorem for spherical 
harmonics and the Jacobi- Anger expansion of a plane wave given 
by (1) and (2) respectively, we find 

^mu.s' = 4;r ^ i^ jt{27:\\u\\) Y^Ju) Y,^(s') . (11) 
Using this result, the integral contained in (10) becomes 

f e-^2^"-^' Y,,n(t) dn(t) = 47T i-iY M27t\\u\\) YUu) , 

where we have noted the orthogonality of the spherical harmonics. 
The harmonic representation of the visibility function then reads 

n^(M) = 4;r ^(-i)' ;V(2;r||M||) Ff„(M) (A' • /%,„ . (12) 

€m 

This expression has been derived independently by Ng (2001) to 
study properties of the angular power spectrum recovered from in- 
terferometric observations of the CMB, however slightly different 
definition are adopted to those used here. 

Computing visibilities using (12) ensures that full- sky con- 
tributions from the beam-modulated intensity function to the vis- 
ibility integral are incorporated, thus allowing for arbitrarily wide 



fields of view and beam sizes. The complexity of evaluating (12) 
for a single baseline u scales as (9(^max^), where 0(£jriax^) ~ 0(N) 
(as discussed in more detail in Sec. 4.1). Low-resolution simula- 
tions of full- sky interferometric observations are computed in this 
manner in Sec. 4. 



3.2.2 Image reconstruction 

In theory one may recover a full- sky synthesised image from the 
spherical harmonic representation of the visibility function by in- 
tegration over a spherical surface in . We describe this approach 
before discussing its limitations. 

It is possible to recover the beam-modulated intensity func- 
tion by taking the spherical harmonic transform of the visibility 
function over a spherical surface at radius By substituting the 
harmonic representation of the visibility function given by (12) into 
the spherical harmonic transform of the visibility one finds 

r ^(«) y;Ju) dQ(u) = An {-{)' jt{2n\\u\\) {A' • (13) 

where we have noted the orthogonality of the spherical harmonics. 
In theory, (13) can be used to recover the harmonic coefficients of 
the beam-modulated intensity, from which the real space function 
can be reconstructed easily. In this framework, one recovers the 
beam-modulated intensity on the entire sky, arbitrarily far from the 
interferometer pointing direction. However, there are a number of 
limitations of this approach that render it infeasible in practice. 

To recover the harmonic coefficients {A} • f)^^ from (13) we 
require that 7V(27r||M||) is non-zero for the particular value of £ 
and the radius \\u\\ considered. Firstly, let us attempt to recover 
(A^ • f)^^^ for all £ and m from a spherical sampling of the visi- 
bility function at a single radius In theory this is possible since 
the zth zeros of the spherical Bessel function are strictly monotoni- 
cally increasing with order £ (Liu & Zou 2007), hence by choosing 
a value of 27t\\u\\ that lies between two zeros of identical zero-order 
z and adjacent function-order ^, we ensure that we avoid all zeros 
of the spherical Bessel functions. Nevertheless, as the order £ devi- 
ates from the adjacent values chosen, the spherical Bessel functions 
will become arbitrarily closer to zero. Consequently, any numerical 
attempts to recover (A^ • f)^^^ using this procedure will be unstable. 
A full u sampling of the visibility function in is therefore re- 
quired to recover the full- sky beam-modulated intensity function. 
For each value of £, sampled spherical surfaces with different radii 
should be used. The value 27t\\u\\ used for a particular £ should be 
chosen to ensure that it does not lie on or near the zeros of the 
spherical Bessel function 7V(2;7r||M||), and ideally at a low-order lo- 
cal extremum of the spherical Bessel function. 

We have demonstrated that it is possible to recover the beam- 
modulated intensity over the entire sky in theory, however to do 
this we require full sampling of the visibility function in R^. For 
interferometric observations, as the field of view rotates over the 
sky with time we sample the visibility function for various values 
of the baseline in local coordinates. Typically these samples lie on 
a series of wv-tracks in the u = (w, v, w) space for values of w close 
to zero. For low pointing directions we recover samples at values of 
w further from zero, but we are unlikely in practice to recover suf- 
ficient sampling of the visibility function to reconstruct the beam- 
modulated intensity function over the entire sky. Consequently, we 
do not advocate the use of the spherical harmonic representations 
derived above for image recovery. 
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3.3 Spherical Haar wavelet space representation 

The main challenges that arise from the full- sky interferometry for- 
mulations in both real and spherical harmonic space are related to 
the localisation of signal characteristics. In real space we achieve 
good spatial localisation but poor frequency localisation; in spher- 
ical harmonic space we achieve good frequency localisation but 
poor spatial localisation. Both the primary beam and the sky in- 
tensity functions are characterised by spatially localised high fre- 
quency content, thus neither real nor harmonic space bases provide 
an efficient representation of these functions. However, these func- 
tions can be represented efficiently in a wavelet basis due to the 
simultaneous spatial and scale localisation aff'orded by a wavelet 
analysis (as discussed in Sec. 2.2). We derive here the SHW rep- 
resentation of the interferometry integral and discuss the implica- 
tions of this representation for computing visibilities and for image 
reconstruction. 

3.3.1 Computing visibilities 

Consider the SHW representation of the beam-modulated intensity 
function 

^/Q-l 7-1 A^;-l 2 

(A^ • = £ Xj,, <pj,,{t) + X Z Z ' (14) 

k=Q j=Jo k=0 m=0 

where Aj^^k and yjj^ are the scaling and wavelet coefficients of the 
beam-modulated intensity function respectively, and the SHW rep- 
resentation of the plane wave 

k=0 j=Jo k=0 m=0 

where rjj^^k and SJj^ are the scaling and wavelet coefficients of the 
plane wave respectively. Substituting these expansions into the lo- 
cal coordinate visibility integral given by (9), we obtain 

"^("^^ Z Tj '^'o,kr]jo,k'(u) I (Pj,,k(s')(/>j^,k<s')dQis') 

k=0 k'=0 

+ Z Z Z Z^^o,^ 0yo,^(^V/;,,(^V^^(^^) 

k=0 f=jQ k'=0 m'=0 
j-lNj-l 2 ^/o-l 

+ ZZZ Z ^r^^^o,^") i//l,(s')cpj,,As')daCs') 

j=jQ k=0 m=0 k'=0 

j=7o ^=0 m=0 j'=jQ k'=0 m'=0 

Noting the orthogonality of the scaling functions and wavelets this 
reduces to 

^/Q-l 7-1 A^;-l 2 

^(u) = £ Aj,,, rjj,,,(u) + X Z Z ^7,k(^^ ■ (1^) 

k=0 j=Jo k=0 m=0 

Applying (16) naively to compute visibilities in the full-sky setting 
is no more efficient than computing visibilities in real or harmonic 
space and scales as 0(N) (since the multiresolution SHW decom- 
position contains exactly as many scaling and wavelet coefficients 
as the number of pixels on the original sphere). Furthermore, to 
compute (16) one must first compute the wavelet coefficients of the 
plane wave though (15), for each baseline u. The spherical harmon- 
ics define a natural harmonic representation on the sphere, arising 



from the solution to Laplace's equation in spherical coordinates. 
Consequently, the expansion of the plane wave e^^^" can be per- 
formed analytically in the spherical harmonic basis through (11). 
Unfortunately this is not possible in the SHW basis and the scaling 
and wavelet coefficients of the plane wave must be computed nu- 
merically, introducing an additional overhead when computing vis- 
ibilities in SHW space. However, these disadvantages are more than 
off'set by the efficient representation of the beam-modulated inten- 
sity function in SHW space. Beam-modulated intensity functions 
are likely to contain localised high-frequency content and hence 
will be extremely sparse in the wavelet basis, i.e. a large number 
of wavelet coefficients will be zero. If we consider only the non- 
zero wavelet coefficients when computing (16), visibilities can be 
computed at substantially lower computational expense. Due to the 
pairing of approximation and detail coefficients in (16), it is also 
only necessary to compute the scaling and wavelet coefficients of 
the plane wave that correspond to the non-zero coefficients of the 
beam-modulated intensity function, thus reducing the overhead dis- 
cussed above. Although it is not possible to determine exactly how 
(16) scales with the resolution of the pixelised spheres considered, 
since this depends on the SHW sparsity properties of the particu- 
lar beam-modulated intensity function considered, the complexity 
of computing (16) is found (see Sec. 4) to something scale like 
^(Cax"), where w < 1. 



In practice, not only are many of the SHW coefficients of the 
beam-modulated intensity identically zero, but many are also very 
close to zero. These wavelet coefficients contain minimal informa- 
tion content and can be set identically to zero without introduc- 
ing substantial errors in the representation of the original beam- 
modulated intensity function on the sphere. For practical imple- 
mentations a strategy is required to determine those wavelet coef- 
ficients that are sufficiently close to zero that they can be safely 
ignored. We develop two such strategies. The first strategy involves 
simple hard thresholding, so that all wavelet coefficients below a 
threshold (in absolute value) are ignored. The threshold value is 
determined by specifying the proportion of wavelet coefficients to 
retain. Typically less than one percent of the wavelet coefficients 
can be kept while still ensuring visibilities are computed accurately, 
as we shall see in Sec. 4. This strategy treats all of the wavelet coef- 
ficients identically. However, wavelet coefficients on coarser levels 
(lower j) are defined over a larger portion of the sky and so con- 
tain more information content. One should therefore favour keeping 
wavelet coefficients at coarser levels over finer levels. Our second 
strategy for determining the wavelet coefficients to retain is based 
on hard thresholding with an annealing strategy to determine the 
threshold value separately for each level j. The proportion of coef- 
ficients to retain at the finest level is specified and this proportion 
is increased quadratically as ones progresses to coarser levels. We 
also experimented with linear and exponential annealing strategies, 
however the quadratic strategy was most successful in characteris- 
ing the importance of wavelet coefficients across levels. 

In Sec. 4 low-resolution simulations of full-sky interferomet- 
ric observations are computed using the SHW method outlined 
here. Due to the superior computational efficiency of the SHW 
method compared to the real and spherical harmonic space ap- 
proaches, it is also possible to perform high-resolution simulations 
of interferometric observations using this method; these are also 
presented in Sec. 4. 
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3.3.2 Image reconstruction 

In the SHW formulation of full-sky interferometry given by (16), 
the scaling and wavelet coefficients of the beam-modulated inten- 
sity function are linearly related to the visibilities by the scaling and 
wavelet coefficients of the plane wave. The plane wave is oscilla- 
tory and does not contain localised high frequency content, hence 
the plane wave is not likely to be sparsely represented in the SHW 
basis. Consequently, the use of (16) for wide field of view image re- 
construction is much more well posed than the spherical harmonic 
representation of the visibility integral. Furthermore, if the primary 
beam is known, which is often the case, then one need only attempt 
to recover the wavelet coefficients of the beam-modulated intensity 
function that correspond to the non-zero coefficients of the primary 
beam, i.e. the beam tells us exactly which wavelet coefficients to 
recover and all others can safely be assumed to be zero. 

SHW image reconstruction on the full- sky therefore involves 
inverting the linear system given by (16). To solve this system is it 
instructive to write it as the matrix equation 





T 








7 



■ A^{u) r , 



where Aj^ and y are vectors of concatenated non-zero scaling and 
wavelet coefficients of the beam-modulated intensity and Jjj^iu) 
and 6(u) are the corresponding coefficients of the plane wave. Com- 
bining the scaling and detail coefficients into a single vector we ob- 
tain r = [Ajq y]^ and \(u) = [^/^(m) S(u)]^. Now assuming that 
we have M visibility observations corresponding to diff'erent base- 
lines, we may write the system of equations that we recover as the 
single matrix equation 



(17) 



where 



and 



a^=[^(M0) ^(Ml) ^(MM)f 

M = [A(mo) A(mi) A(mm)] . 

The overdetermined system (17) may solved in the least squares 
sense for an estimate of the non-zero scaling and wavelet coeffi- 
cients of the beam-modulated intensity function: 



(18) 



One would expect that recovering the beam-modulated inten- 
sity function on the full- sky in this manner would be well posed, 
however to really ascertain the eff'ectiveness of the method outlined 
here it should be implemented and tested. The focus of the current 
article is predominantly on the forward interferometric wide field 
imaging problem. The use of SHWs for wide field of view image 
reconstruction is an interesting application in its own right and we 
intend to develop these ideas further and present the results of var- 
ious experiments in a future work. 



3.4 Incorporating horizon occlusion and variable beams 

In the preceding formulations we have assumed that the full- sky 
is visible. Obviously this is not the case for Earth based interfer- 
ometers that may observe only the hemisphere above the horizon. 
Nevertheless, if the interferometer pointing direction is relatively 
close to the North pole of the observable hemisphere and the beam 
is sufficiently small that it is zero in the southern hemisphere, then 



the full-sky formulation is appropriate. Furthermore, we have as- 
sumed that the primary beam is fully defined in local coordinates 
and is independent of the interferometer pointing direction on the 
sky. This is unlikely to be the case for real aperture array interfer- 
ometers. In this section we formulate full- sky interferometry in the 
setting where the beam is sufficiently large that it may be non-zero 
below the horizon and we discuss extensions to incorporate primary 
beams that also depend on pointing direction. 

In order to incorporate horizon occlusion and variable beams 
we must introduce a third coordinate system. In addition to the lo- 
cal coordinate system defined relative to the interferometer point- 
ing direction and the global coordinate system of the celestial sky, 
it is necessary to introduce an Earth-based coordinate system. We 
define the Earth-based coordinate frame by the right handed set of 
unit vectors {^1,^2,^3)- The Earth-based coordinate system is re- 
lated to the celestial sky coordinate system by a time varying ro- 
tation % that relates the orientation of the celestial sky relative to 
the Earth. Let denote the corresponding time varying 3 x 3 ro- 
tation matrix relating these coordinate systems. The unit vectors 
defining these coordinate systems are then related by hi = Rt^i. We 
adopt the convention that these coordinate frames are aligned for 
time ? = 0, i.e. R^=o = h where I is the identity matrix. A vector in 
Earth-based coordinates is related to a vector in celestial sky coor- 
dinates by = R~^5", where the superscript e denotes Earth-based 
coordinates. In particular, the interferometer pointing direction in 
Earth-based coordinates traces a fixed point in the celestial sky Sq 
over time: sl(t) - R^^^O' pointing direction in Earth-based 

coordinates is time dependent. The local coordinate system is now 
defined relative to 5o(0 and is related to the Earth-based coordinate 
system by the rotation = 9^((p^, ^^,0), with corresponding rota- 
tion matrix Rq, where (6^, ip^) are the spherical coordinates of sl(t). 
For notational simplicity the time dependence of "Rq and Rq is left 
implicit. 

Horizon occlusion may be modelled by incorporating a binary 
horizon function that is unity above the Earth's horizon and zero 
below. This horizon function is represented most naturally in the 
Earth-based coordinate system H^(s^) e L^(S^, dQ) and is defined 
by 



1 if r • ^3 > 

otherwise . 



Representing each function in its natural coordinate system the in- 
terferometer visibility integral may be written 



JS2 



(19) 



where the binary horizon function is included to exclude contribu- 
tions to the visibility integral from below the horizon. It is again 
convenient to represent all functions in the visibility integral in a 
consistent coordinate system. In local coordinates the horizon func- 
tion becomes 

i.e. -Kq^W, and the source intensity function becomes 



r(r) = r(R,Ro5^) = (K^'n;'p){s') = i\s') , 

i.e. f = Kq^K'^F. Any azimuthal 7 rotation of % will render P 
time dependent due to the fixed azimuthal rotational component of 

The local coordinate version of the visibility integral then 
follows trivially from (19). The techniques outlined in Sec. 3.1 
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through Sec. 3.3 can then be appUed directly by replacing the beam- 
modulated intensity function (A^ • with the horizon-beam- 
modulated intensity function (A^ • ■ = A\f)H\sY(s^). 
The time dependence of and f now precludes any full- sky 
reconstruction of the beam-horizon-modulated intensity function 
even in theory, let alone in practice. The interferometer baseline 
in local coordinates is given by m = Rq^m^, where is the natural 
representation of the baseline in Earth-based coordinates. Note that 
u is time dependent due to the time dependence of Rq. It was not 
possible to define u explicitly previously since we did not have a 
natural coordinate system to represent the interferometer baseline 
(due to the absence of the Earth-based coordinate system). 

If the primary beam function depends on the pointing direction 
in Earth-based coordinates then it also becomes time dependent, i.e. 
A^{f, slit)). In this case the representations of the visibility integral 
derived previously in real, spherical harmonic and SHW space all 
still hold, however A^(f) must be recomputed for each pointing di- 
rection. In this framework it is also possible to include more com- 
plicated horizon eff'ects in the primary beam, rather than the simple 
masking adopted by including the binary horizon function. A sim- 
ulated interferometer visibility observation incorporating horizon 
occlusion and a variable beam would be performed as follows: 

(i) compute the Earth-based pointing direction by sl(t) = R^^^O' 

(ii) compute the interferometer baseline in local coordinates by 
u = R-^M^ 

(iii) compute the intensity function in local coordinates by 

(iv) compute the horizon function in local coordinates by 

(v) compute the primary beam function for the given pointing 
direction by A^it, sl{t))\ 

(vi) compute the observed visibility "Viu) by the method of 
choice using either (9), (12) or (16); 

(vii) repeat steps (i)-(vi) for each observation time t. 

Although the beam may be sufficiently large that a small patch 
approximation is not valid, it may often be the case that it is not so 
large that horizon occlusion must be included when computing vis- 
ibilities. Consequently, although we have presented the most gen- 
eral formulation of full- sky interferometry in this section by care- 
fully modelling horizon occlusion and including variable beams, it 
may often be appropriate to neglect these effects and simply follow 
the formulation outlined previously. 



4 SIMULATIONS 

We perform simulations of interferometric observations, including 
full- sky contributions, in order to validate the full- sky interferome- 
try formulations presented in Sec. 3. Firstly, we compare all of the 
methods on low-resolution simulations. Performing full- sky inter- 
ferometric simulations using either the real or spherical harmonic 
space representation is a substantial computational challenge and 
is not feasible for high-resolution simulations. The SHW meth- 
ods ease this computational burden considerably, rendering higher 
resolution simulations feasible. Using SHW methods we also per- 
form high-resolution simulations of full- sky interferometric obser- 
vations. Before presenting these low- and high-resolution simula- 
tions we begin with a brief discussion of the practical considera- 
tions that must be taken into account. 



4.1 Practical considerations 

In this article we focus predominantly on the forward wide field 
imaging problem that involves simulating the visibilities observed 
by an interferometer when including full-sky contributions. In 
terms of the inverse problem of image reconstruction on wide fields 
of view, we have seen that this problem is not well posed in the 
spherical harmonic representation but is likely to be well posed 
in the SHW space representation. We have outlined a preliminary 
method to perform image reconstruction using SHWs in Sec. 3.3.2. 
In a future work we intend to develop these ideas further and to ex- 
amine the performance of these methods on simulations. However, 
we restrict our attention in the simulations performed in this arti- 
cle to the forward problem of simulating visibilities. Consequently, 
although we consider full- sky effects here when computing visibil- 
ities, we adopt the usual Fourier transform approach for recovering 
images on the tangent plane at the interferometer pointing direc- 
tion. We next relate the parameters of our visibility simulations, 
including the visibilities that must be computed, to the resolution 
and size of the recovered tangent plane image. 

In practice we assume that all signals on the sphere are band- 
limited at ^max- Although the exact number of samples on the sphere 
N required to represent a band- limited signal depends on the pixeli- 
sation of the sphere, for the HEALPix scheme and in general it is of 
order (9(A^) 0{^mJ) (e.g. Driscoll & Healy 1994; Gorski et al. 
2005). The maximum interferometer baseline distance ||M||max for 
which the visibility may be computed accurately is related to the 
harmonic band- limit ^max- For small fields, the approximate rela- 
tionship ^max - 2:7r||M||inax holds (the exact relationship is 2||M||max = 
cot(;7r/^niax); Hobson & Magueijo 1996). In order to ensure Nyquist 
sampling is satisfied when reconstructing images on small patches, 
for square images one requires that /^max = V^i"^age/(2wmax), where 
/7max is the width of the image plane, w^ax = ll«llmax/ V2 and Mmage 
is the number of pixels of the reconstructed image (these results 
generalise to rectangular images trivially). 

In order to image these small patches accurately, a high band- 
limit ^max must be considered. This constraint increases the compu- 
tational requirements of computing visibilities in harmonic space 
significantly. For example, to image an area on the sky of one 
square degree (pmax = 1°) for a 50x50 pixel image (Mmage = 
2, 500), we must consider all harmonic coefficients up to a band- 
limit of ^max - 13, 000. It is then necessary to repeat this computa- 
tion for numerous local coordinate representations of the interfer- 
ometer baseline u as the source traverses the sky. Although an exact 
sampling theorem does not exist on a HEALPix sphere, typically the 
harmonic band- limit of a function sampled on a HEALPix sphere 
is related to the resolution of the sphere through ^^ax = cA^side, 
where c is typically two or three. In order to accurately image these 
small patches in real space while including full- sky contributions 
a high resolution pixelisation of the sphere is therefore required. 
For the imaging example considered above, one must consider all 
pixel values of a HEALPix sphere at resolution A^side ^ 512. As dis- 
cussed previously, the advantage of the SHW representation is its 
simultaneous localisation of signal content in scale and position. 
Computing full- sky visibilities using the SHW method adapts to 
the sparse representation of the beam-modulated intensity function 
in the SHW basis, thus providing substantial computational savings 
over the real and spherical harmonic space methods. We next ap- 
ply all of these methods to low-resolution simulations to compare 
their performance. Finally, let us note that the visibility computa- 
tion for each baseline u is independent and implementations of all 
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of these methods may be paralleUsed easily to spread the computa- 
tional load. 



4.2 Low-resolution simulated synchrotron observations 

In this section we simulate low-resolution observations of syn- 
chrotron emission made by an idealised interferometer. We sim- 
ulate visibilities using all of the the full- sky interferometry frame- 
works outlined in Sec. 3, ensuring that contributions due to a large 
primary beam are included. For simplicity, we do not incorporate 
in these simulations the extensions described in Sec. 3.4 that al- 
low horizon occlusion and variable beams. The motivation of these 
simulations is to demonstrate and validate the application of the 
full- sky interferometry formulations presented previously and to 
compare the performance of the methods. Our implementation of 
these methods have yet to be parallelised and all timing tests pre- 
sented here and subsequently are performed on a laptop with a sin- 
gle 2.2GHz Intel Core 2 Duo processor and 2GB of memory. 

For the background intensity map we use the full- sky syn- 
chrotron template recovered from the 3 -year Wilkinson Microwave 
Anisotropy Probe (WMAP) observations (Hinshaw et al. 2007). A 
detailed description of the construction of this template is given 
by Hinshaw et al. (2007), however for our purpose we merely con- 
sider this as the full- sky background intensity map F{s^) to be ob- 
served by our idealised interferometer. This synchrotron map is 
available for download from the Legacy Archive for Microwave 
Background Data Analysis (LAMBDA).^ In order to remove high- 
frequency content of this map, we smooth the data with a Gaussian 
kernel with full-width-half-maximum of FWHMs = 1.7° (since we 
only perform low-resolution simulated observations of these data). 
The resulting smoothed, full- sky synchrotron map is illustrated in 
Fig. 5. We simulate observations of the extended source located at 
^0 = (^S'*^o) = (84.0°, 76.5°). In order to compute full-sky visi- 
bility contributions, it is necessary to convert the global intensity 
function to its local coordinate version centred on the interferome- 
ter pointing direction. This is performed through a rotation by 
as outlined in Sec. 3.1 . The local coordinate intensity function I^is^) 
is illustrated in Fig. 6 (a). The extended source to be observed is 
clearly visible at the north pole of this map. The interferometer 
beam used in these simulations is given by a wide Gaussian with 
FWHMb - 18° and is illustrated in Fig. 6 (b). 

Given the mock data and configuration discussed previously, 
we simulate visibilities through direct quadrature on the sphere 
given by (9), the spherical harmonic space representation given by 
(12) and the SHW space representation given by (16). We assume 
complete uv coverage and consider a baseline limit of Wmax = 30, 
corresponding to ^max - 270, and reconstruct an image of size 
Mmage = 20 X 20. The smoothed synchrotron map is represented 
by a HEALPix sampled sphere at resolution A^side = 256 and so is 
sampled sufficiently for these simulations. Once the visibilities are 
simulated, a synthesised image is reconstructed simply by taking 
the inverse Fourier transform. Nyquist sampling dictates that the 
synthesised image corresponds to a ~ 20° square patch. This im- 
age reconstruction relies on a tangent plane approximation which 
will not be valid for the large field of view considered. Neverthe- 
less, this simple approach to reconstruction is sufficient to demon- 
strate the validity of our full- sky interferometry framework in this 
restricted low-resolution setting. The application of the standard 
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(a) Mollweide projection 




(b) Globe 

Figure 5. Full-sky synchrotron map observed by WMAP and smoothed 
with a Gaussian kernel of FWHMs = 1.7°. This synchrotron map pro- 
vides the full-sky background intensity map for our low-resolution simu- 
lated observations and is shown here in the global coordinate frame defined 
by Galactic coordinates. 



Fourier transform here to synthesise an image from the simulated 
visibilities is performed for visual verification only. 

The direct projection of the full- sky beam-modulated inten- 
sity function onto the tangent plane at the interferometer pointing 
direction is illustrated in Fig. 7 (a). The images reconstructed from 
full- sky visibilities computed using all of the methods outlined in 
Sec. 3 are illustrated in the remaining panels of Fig. 7 (upsam- 
pled to a 60 X 60 pixel image for visualisation). One would ex- 
pect the tangent plane image of Fig. 7 (a) to diff'er slightly from 
the reconstructed images since full- sky contributions due to the 
large beam are incorporated when simulating visibilities, however 
a flat-patch approximation is assumed when synthesising the im- 
age using the standard Fourier approach. Furthermore, the images 
synthesised from simulated visibilities contain less high-frequency 
content as a consequence of the band-limit imposed by the interfer- 
ometer baseline limit. Nevertheless, all images are in close agree- 
ment, demonstrating and validating the use of the methods derived 
in Sec. 3 to simulate visibilities observed by an interferometer in 
the full-sky setting. Furthermore, the images contained in panels 
(b), (c) and (d) of Fig. 7, corresponding respectively to the real, 
spherical harmonic and naive SHW methods of computing full- sky 
visibilities are identical (to numerical precision). These images are 
computed using independent implementations, thereby providing 
an additional verification of these methods and implementations. In 
order to realise the advantages of the efficient representation of the 
SHW basis, it is necessary to ignore near-zero valued wavelet co- 
efficients when computing (16). This is not done in the naive SHW 
method but strategies to disregard unimportant wavelet coefficients 
were described in Sec. 3.3.2: constant and annealing based thresh- 
olding strategies were proposed. Reconstructed images based on 
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are already considerably faster at this low-resolution, thus render- 
ing high-resolution simulations feasible. 




(a) Synchrotron map 




(b) Gaussian beam profile with FWHMb - 18° 

Figure 6. Full- sky synchrotron and beam maps in local coordinates. We 
consider the source at = (6q,(Pq) = (84.0°, 76.5°). The synchrotron map 
in global coordinates is rotated hy^R^^ = 'RiO, -6^, -(p^) to convert to local 
coordinates. Full globes are shown on the left; zoomed images of globes 
about the north pole are shown on the right. Note that the beam profile takes 
the value of unity in the centre and -0.3 at the boundary of the zoomed 
image. 



these methods are illustrated in panels (e) and (f) of Fig. 7. These 
approaches do introduce a small error compared to the images of 
panels (b) through (d) of Fig. 7, but the controlled introduction of 
these errors provide substantial computational savings. 

In Table 1 we compare the performance of the methods used 
to compute visibilities in the simulations presented in Fig. 7. The 
complexity of each method was discussed in Sec. 3 and we collate 
these complexities here and relate them to the maximum baseline 
of the interferometer ||w||max- For the exact full-sky interferometry 
formalisms all coefficients are used, be they pixel values, spherical 
harmonic coefficients or scaling and wavelet coefficients. However, 
the thresholded SHW strategies require only a small proportion of 
wavelet coefficients due to the sparse representation of the beam- 
modulated intensity function in the SHW basis. Typically, less than 
one percent of the wavelet coefficients can be retained while intro- 
ducing only minimal errors. It is apparent that the annealed thresh- 
olding strategy is slightly superior to the constant thresholding ap- 
proach, as one would expect since the annealing approach is more 
likely to retain the detail coefficients of coarser levels that inher- 
ently contain more information content. The execution time tests 
presented in Table 1 show that the thresholded SHW methods are 
considerably faster than the other methods at this resolution. How- 
ever, the true advantages of the thresholded SHW methods will only 
become apparent at higher resolution simulations. At higher resolu- 
tions, the real, spherical harmonic and naive SHW methods all scale 
like (9(||M||inax^)- The already slow performance of these techniques 
and their poor scaling render these methods computationally infea- 
sible for high-resolution problems. The thresholded SHW methods 
have much better scaling properties (as discussed in Sec. 3.3.1) and 



4.3 High-resolution simulated dust observations 

In this section we simulate high-resolution observations of Galactic 
dust emission made by an idealised interferometer. Full-sky visi- 
bilities are simulated using the annealing based thresholded SHW 
method only since computational considerations preclude the ap- 
plication of real and spherical harmonic representations, and the 
annealing based SHW method was demonstrated in Sec. 4.2 to be 
superior to the constant thresholding method. We subsequently re- 
fer to the annealing based SHW method for simulating full- sky vis- 
ibilities as the fast SHW method. 

The synchrotron map considered previously does not contain 
enough high-frequency content for the high-resolution simulations 
performed here. For the background intensity map we therefore use 
the 94GHz FDS map of predicted submillimeter and microwave 
emission of diffuse interstellar Galactic dust (Finkbeiner et al. 
1999). This predicted map is based on the merged Infrared Astron- 
omy Satellite (IRAS) and Cosmic Background Explorer Diffuse In- 
frared Background Experiment (COBE-DIRBE) observations pro- 
duced by Schlegel et al. (1998). An undersampled version of the 
FDS map is available from LAMBDA as a HEALPix sphere at res- 
olution A^side = 512. The FDS map that we assume as our full-sky 
background intensity map P(s^), to be observed by our idealised 
interferometer, is illustrated in Fig. 8 (a) and (b) in global coordi- 
nates. We simulate observations of the extended source located at 
Sq = (Oq,(Pq) = (108.0°, 0.0°). It is necessary to convert the global 
intensity function to its local coordinate version through a rotation 
by in order to compute visibilities. The local coordinate inten- 
sity function is illustrated in Fig. 8 (c) and (d). The extended 
source to be observed is clearly visible at the north pole of this 
map. The interferometer beam used in these simulations is given 
by a Gaussian with FWHMb ^ 2.9°. 

Given the mock data and configuration discussed above, we 
simulate visibilities using the fast SHW method (i.e. the annealing 
based thresholded SHW approach). We assume complete uv cover- 
age and consider a baseline limit of Wmax = 100, and reconstruct an 
image of size Mmage = 20 x 20. Following the approach performed 
in Sec. 4.2, a synthesised image is reconstructed from simulated 
visibilities simply by taking the inverse Fourier transform. The lim- 
itations of this approach were discussed previously, however it is 
suitable for visual verification of the simulations performed in this 
article. Nyquist sampling dictates that the synthesised image corre- 
sponds to a ~ 5.7° square patch. 

The direct projection of the full- sky beam-modulated intensity 
function onto the tangent plane at the interferometer pointing direc- 
tion is illustrated in Fig. 9 (a). The beam-modulated intensity image 
reconstructed from the full-sky simulated visibilities is illustrated 
in Fig. 9 (b) (upsampled to a 60x60 pixel image for visualisation). 
In this simulation visibilities are computed using 0.023 percent 
of the wavelet coefficients of the beam-modulated intensity func- 
tion only. Again, one expects the tangent plane and reconstructed 
images to differ slightly since full- sky contributions are incorpo- 
rated when simulating visibilities, however a flat-patch approxima- 
tion is assumed when synthesising the image using the standard 
Fourier approach. Furthermore, the fast SHW method itself intro- 
duces some small error when discarding those wavelet coefficients 
with minimal information content and the interferometer baseline 
limit also restricts the high-frequency content of the reconstructed 
image. Nevertheless, the tangent plane and reconstructed image are 
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Figure 7. Beam-modulated synchrotron intensity images for a -20° square patch. The image shown in panel (a) is constructed by projecting the full-sky 
beam-modulated intensity onto the tangent plane at the interferometer pointing direction defined by the coordinate system p = (p, q) (as defined in Sec. 3.1). 
The images shown in the remaining panels are constructed by simulating visibilities in the full-sky setting using all of the methods described in Sec. 3, 
followed by a standard inverse Fourier transform to recover the synthesised image. The image in panel (a) is expected to diff'er to the other images since 
full-sky contributions due to the large beam are incorporated when simulating visibilities, however a flat-patch approximation is assumed when synthesising 
images using the standard Fourier approach. Nevertheless, all images are relatively similar, demonstrating and validating the use of the visibility formulations 
derived in Sec. 3 to simulate visibilities observed by an interferometer in the full-sky setting. The reconstructed images shown in panels (b), (c) and (d) are 
identical (to numerical precision) as expected, while the reconstructed images shown in panels (e) and (f) difl'er to these very slightly since some information 
is discarded when thresholding the wavelet coefficients in order to increase the speed of computations. 



Table 1. Comparison of the performance of the methods derived in Sec. 3 for simulating interferometric observations in the full-sky setting. The results 
presented in this table correspond to the simulations illustrated in Fig. 7. Execution time tests were performed on a laptop with a 2.2GHz Intel Core 2 Duo 
processor and 2GB of memory. 



Method 


Complexity 0(||M||max") 


Coefficients retained^ 


Execution time 


Direct quadrature 


n = 2 


100.00% 


207.6s 


Spherical harmonic 


n = 2 


100.00% 


263.7s 


Naive SHW 


n = 2 


100.00% 


238.9s 


Thresholded SHW (constant threshold) 


n< 1 


0.70% 


75.8s 


Thresholded SHW (annealing strategy) 


n<l 


0.35% 


73.0s 



^ Note that the number of coefficients retained is only relevant for the thresholded SHW methods; all other methods require all coefficients (be they pixel 
values, spherical harmonic coefficients or SHW coefficients). 



in close agreement, demonstrating and validating the application of 
the fast SHW approach to simulate full- sky visibilities observed by 
an interferometer in a high-resolution setting. 

The high-resolution simulations performed here using the fast 
SHW method required an execution time of 289.6s. Based on the 
scaling relationships of the real and spherical harmonic space meth- 
ods, and their execution times on the low-resolution simulations 



performed in Sec. 4.2, we estimate the execution time of these 
methods when applied to these high-resolution simulations to be 
~ 3000s. As the resolution of the problem (i.e. the baseline limit 
Umax) and the number of visibilities to be computed (i.e. Mmage 
in this notation) increase, this difference will become even more 
marked. The fast SHW method is therefore essential to compute re- 
alistic simulations of visibilities when incorporating full- sky con- 
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(a) MoUweide projection in global coordinates 




(b) Globe in global coordinates 




(c) Globe in local coordinates 




(d) Zoomed globe in local coordinates 



Figure 8. Full-sky 94GHz FDS map of predicted submillimeter and mi- 
crowave emission of diffuse interstellar Galactic dust. This map provides 
the full-sky background intensity map for our high-resolution simulated ob- 
servations. Panels (a) and (b) show the full-sky map in the global coordi- 
nate frame defined by Galactic coordinates. We simulate observations of the 
source at = (6^, (p^ = (108.0°, 0.0°). A rotation by "R'^ = ^(0, -6^, -(p^ 
is performed to convert the map to the local coordinate versions illustrated 
in panels (c) and (d), which are centred on the north pole. 




(b) Simulated full-sky interferometric image 



Figure 9. Beam-modulated Galactic dust intensity images for a -5.7° 
square patch. The image shown in panel (a) is constructed by projecting 
the full- sky beam-modulated intensity onto the tangent plane at the inter- 
ferometer pointing direction defined by the coordinate system p = (p, q) 
(as defined in Sec. 3.1). The image shown in panel (b) is constructed by 
simulating visibilities in the full-sky setting using the fast SHW method, 
followed by a standard inverse Fourier transform to recover the synthesised 
image. The images in panels (a) and (b) are expected to diff'er since full- sky 
contributions due to the large beam are incorporated when simulating visi- 
bilities, however a flat-patch approximation is assumed when synthesising 
the image using the standard Fourier approach. Furthermore, the fast SHW 
method itself introduces some small error when discarding those wavelet 
coefficients with minimal information content. Nevertheless, the two im- 
ages are relatively similar, demonstrating and validating the use of the fast 
SHW method to simulate full-sky visibilities in a high-resolution setting. 

tributions. In future, we intend to parallelise our implementation 
of the fast SHW method and apply it to produce realistic simula- 
tions of full-sky interferometer observations, including incomplete 
uv coverage and other more realistic assumptions. 



5 CONCLUSIONS 

Next generation interferometers will have very large fields of view, 
which poses a number of imaging challenges. The usual interfer- 
ometric approach to simulate visibilities or reconstruct images is 
based on standard Fourier imaging, which relies on a tangent plane 
approximation that is valid only for a small field of view. This ap- 
proach is inappropriate for the wide fields of view proposed for 
the next generation of interferometers. In this article we have for- 
mulated interferometry in the full- sky setting, incorporating con- 
tributions over the entire sky to ensure that contamination due to 
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wide sidelobes of primary beams is not neglected. Full- sky inter- 
ferometry formulations have been developed in real, spherical har- 
monic and SHW spaces. The SHW formulation proves the most 
advantageous approach due to the efficient representation of the 
spatially localised high-frequency content typical of primary beams 
and sky intensity functions in the wavelet basis. A corresponding 
fast SHW method was developed to simulate full- sky visibilities. 
The fast SHW method exploits the sparsity of the wavelet repre- 
sentation by discarding all wavelet coefficients that contain mini- 
mal information content. A quadratic annealing based thresholding 
strategy was developed to determine those wavelet coefficients that 
may be safely discarded. The resulting fast SHW method has su- 
perior computational scaling properties and may be performed at a 
substantially lower computational cost than the real and spherical 
harmonic space methods for simulating visibilities when including 
full- sky contributions. 

The primary focus of this article is the use of these full- sky 
interferometry representations to simulate visibilities, however we 
also briefly discussed implications for the reconstruction of images 
on wide fields of view. Although it is possible in theory to recon- 
struct full- sky images in the spherical harmonic representation, this 
requires full sampling of the visibility function in 'M? and hence is 
not likely to be well posed in practice. However, image reconstruc- 
tion on wide fields of view is likely to be well posed in SHW space 
and we outlined a preliminary method to perform wide field image 
reconstruction. The use of SHWs for image reconstruction is an in- 
teresting application in its own right and we intend to develop these 
ideas further and to evaluate their eff'ectiveness in a future work. 

To demonstrate the application of our full- sky interferometry 
representations for simulating visibilities, and to compare the vari- 
ous methods, we simulated low-resolution full-sky interferometric 
observations of synchrotron emission using all of the methods. The 
real, spherical harmonic and naive SHW methods recovered iden- 
tical images of the extended source observed by our idealised in- 
terferometer, where the naive SHW method does not exploit sparse 
representations in the wavelet basis. The fast SHW method does 
exploit sparse representations, introducing a small error when dis- 
carding those wavelet coefficients that contain minimal information 
content, but consequently the speed of computations in increased 
substantially. Typically less than one percent of wavelet coefficients 
are retained in these low-resolution simulations, reducing the exe- 
cution time of simulations by a factor of approximately three com- 
pared to the other methods (while introducing only minimal errors 
in reconstructed images). However, the true advantages of the fast 
SHW method only become apparent on high-resolution simulations 
due to its superior scaling properties. 

The already slow performance of the real and spherical har- 
monic space approaches to simulating full-sky visibilities, and 
their poor scaling properties, render these methods computation- 
ally infeasible on higher resolution problems. Computing full- sky 
visibilities using the fast SHW method adapts to the extremely 
sparse representation of the beam and sky intensity function in 
the wavelet basis, thus easing the computational burden of simu- 
lations to the extent that they are rendered computationally feasi- 
ble at high-resolutions. Using this method high-resolution interfer- 
ometric observations of diff'use interstellar Galactic dust were simu- 
lated, demonstrating and validating the use of the fast SHW method 
for simulating visibilities in a high-resolution setting. In this exam- 
ple only 0.023 percent of wavelet coefficients were retained and 
the execution time of simulations was estimated to be greater than 
an order of magnitude faster than simulations based on the real or 
spherical harmonic space full- sky interferometry formulations. 



Now that it is possible to simulate the visibilities observed 
by an interferometer when including contributions over the entire 
sky, a number of related studies may be performed. Firstly, we in- 
tend to parallelise our implementations and develop more realistic 
simulations of full-sky interferometer observations, including in- 
complete uv coverage and other more realistic assumptions. Using 
these simulations, we then intend to study the eff'ect of contamina- 
tion in realistic interferometric observations due to wide sidelobes 
of primary beams. Secondly, our full- sky interferometry formula- 
tion allows one to simulate visibilities for all values ofu - (u, v, w), 
and not only values of u where w = 0, as is the case with the stan- 
dard Fourier transform approach. Such simulations allow one to 
evaluate the performance of the faceting (Cornwell & Perley 1992; 
Greisen 2002) and w-projection (Cornwell et al. 2005) approaches 
to wide field image reconstruction on simulations where the ground 
truth is known. Thirdly, we have highlighted an interesting alter- 
native approach to wide field image reconstruction based on the 
SHW representation that warrants further study. We intend to pur- 
sue all three of these applications in future work. The ability to 
perform realistic high-resolution simulations of interferometric ob- 
servations that include full- sky contributions, afforded by our fast 
SHW formulation, allows many new studies to be performed that 
will be important to the design of next generation interferometers 
and to the development of algorithms to analyse the data generated 
by these instruments. 
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